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Abstract 

We explore a two-dimensional kinematic solar dynamo model in a full sphere, based on 
the helioseismically determined solar rotation profile and with an a effect concentrated 
near the solar surface, which captures the Babcock-Leighton idea that the poloidal field 
is created from the decay of tilted bipolar active regions. The meridional circulation, 
assumed to penetrate slightly below the tachocline, plays an important role. Some doubts 
have recently been raised regarding the ability of such a model to reproduce solar-like 
dipolar parity. We specifically address the parity issue and show that the dipolar mode is 
preferred when certain reasonable conditions are satisfied, the most important condition 
being the requirement that the poloidal field should diffuse efficiently to get coupled 
across the equator. Our model is shown to reproduce various aspects of observational 
data, including the phase relation between sunspots and the weak, diffuse field. 

1 Introduction 

Ever since Parker (1955) formulated the turbulent dynamo theory, varieties of very dif- 
ferent dynamo models for the Sun have appeared in the literature. Within the last few 
years, however, a concensus view seems to be emerging as to what should be the basic 



characteristics of a solar dynamo model. The history of how this concensus view emerged 
is discussed in the Introduction of Nandy & Choudhuri (2001), with citations to impor- 
tant original papers. A more detailed account of the history can be found in the review 
by Choudhuri (2003a). We, therefore, begin here by summarizing the main aspects of this 
concensus view without getting into the details of history again. 

One important ingredient of a solar dynamo model is the differential rotation, which 
has now been mapped by helioseismology. The toroidal magnetic field must be produced 
by the stretching of poloidal field lines primarily within the tachocline — the region of 
concentrated vertical shear at the base of the solar convection zone (SCZ). The toroidal 
field produced in the tachocline would then rise from there due to magnetic buoyancy 
to produce active regions. Simulations of flux rise through the SCZ suggested that the 
toroidal field at the bottom must be of order 10 5 G (Choudhuri & Gilman 1987; Choudhuri 
1989; D'Silva & Choudhuri 1993; Fan et al. 1993; D'Silva & Howard 1993; Caligari et al. 
1995). Since such a strong field is expected to quench the usual mean field a effect (Parker 
1955; Steenbeck et al. 1966), one possibility being considered by many researchers is that 
the poloidal field is generated at the surface from the decay of tilted bipolar region — an 
idea that goes back to Babcock (1961) and Leighton (1969). The poloidal component 
generated at the surface is first advected poleward by a meridional circulation (Wang et 
al. 1989a, 1989b; Dikpati & Choudhuri 1994, 1995; Choudhuri & Dikpati 1999). Finally, 
the poloidal component sinks with the downward flow at the poles and is brought to the 
tachocline where it can be stretched by the differential rotation to generate the toroidal 
field — thus completing the full cycle. Since advection by the meridional circulation plays 
such a crucial role in such a model, we would refer to this model as the circulation- 
dominated solar dynamo or CDSD model. 

Although the above view of the solar dynamo arose by assimilating the ideas of many 
researchers over the years, Wang et al. (1991), Choudhuri et al. (1995) and Durney (1995, 
1996, 1997) were amongst the first to demonstrate the crucial role which meridional 
circulation is expected to play in modern solar dynamo models. While Choudhuri et 
al. (1995) modeled the generation of poloidal field from the decay of active regions by 
introducing a phenomenological a parameter concentrated at the surface, Durney (1995, 
1996, 1997) followed Leighton (1969) more closely to capture this effect by taking two 
flux rings of opposite polarity at the surface. Afterwards, Nandy & Choudhuri (2001) 
have demonstrated that these two approaches give qualitatively similar results. We shall 
develop our models by prescribing an a coefficient concentrated near the surface. Dikpati 
& Charbonneau (1999) and Kiiker et al. (2001) presented CDSD models with solar-like 
internal rotation and with such a specification of a effect. Since dVt/dr has a larger 
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amplitude at the high latitudes within the tachocline (where it is negative) rather than 
at low latitudes (where it is positive), solar-like rotation tends to produce strong toroidal 
fields at high latitudes, in apparent contradiction to the observed fact that sunspots always 
appear at low latitudes. Commenting on this problem, Dikpati & Charbonneau (1999) 
write: "this is an unavoidable inductive effect . . . and no change in model parameters can 
do away entirely with this feature". Kiiker et al. (2001) point out: "All recent dynamo 
models with the observed rotation law are faced with this problem, even when the a effect 
has been strongly reduced in the polar region by the relation a oc sin 2 6 cos 9, as we also 
did". 

Nandy & Choudhuri (2002) have recently shown in a brief communication that this 
problem can be solved by postulating a meridional flow penetrating somewhat deeper 
than hitherto believed. If the meridional flow goes below the tachocline near the poles, 
then the strong toroidal field produced within the tachocline at high latitudes is immedi- 
ately pushed underneath into the convectively stable layers and cannot emerge at the high 
latitudes. The meridional flow then carries this toroidal field through the stable layers to 
low latitudes. There the meridional flow rises and the toroidal flux enters the SCZ to be- 
come buoyantly unstable and produce active regions at low latitudes. Such a penetrating 
meridional flow produces theoretical butterfly diagrams in remarkable agreement with the 
observed butterfly diagram. The conventional wisdom was that the toroidal field which 
forms sunspots at low latitudes must have been produced at the low latitude. With the 
helioseismically determined rotation profile, this seems unlikely and it may well be that 
the toroidal field is actually generated at the high latitude, even though it is not allowed 
to erupt there. 

Recently it has been pointed out by Dikpati & Gilman (2001) and Bonanno et al. 
(2002) that the CDSD model with the Babcock-Leighton mechanism for producing the 
poloidal field near the surface may not give the magnetic configuration with the observed 
parity. Hale's polarity law of bipolar sunspots suggests that the toroidal magnetic field 
is anti-symmetric across the solar equator, implying a dipolar parity. A majority of the 
CDSD models were solved within one hemisphere and the boundary conditions at the 
equator were taken such that the dipolar mode was forced. Dikpati and Gilman (2001) 
solved the dynamo problem in the full sphere and found that the CDSD model with the a 
effect concentrated near the top of the SCZ preferentially excites the quadrupolar mode 
in which the toroidal field is symmetric across the equator — opposite to what is observed. 
Only when the a effect was concentrated near the bottom of the SCZ, they found the 
dipolar parity in conformity with observations. Bonanno et al. (2002) confirmed these 
findings. 
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The aim of the present paper is to provide the technical details of our dynamo models 
not given in the earlier brief paper of Nandy & Choudhuri (2002), to check the parity of 
these dynamo models by extending our code from a hemisphere to a full sphere and to 
address some related issues. In the dipolar mode, the poloidal field lines have to connect 
across the equator. It is necessary for the poloidal field to diffuse efficiently for this to be 
possible. On the other hand, the diffusion of the toroidal field has to be suppressed if we 
want to ensure that the toroidal field has opposite signs on the two sides of the equator. 
Since the toroidal component is much stronger than the poloidal component, we expect 
the turbulent diffusion to be much less effective on the toroidal component than on the 
poloidal component. On using a high diffusivity for the poloidal field and a low diffusivity 
for the toroidal field, we find that the dipolar parity is preferred. There is no need to 
include an additional a effect at the bottom of SCZ to ensure dipolar parity. 

Since the nature of the meridional circulation plays such a crucial role in our model, 
let us make a few comments on it. Within the last few years, helioseismic techniques 
have given us some information about the sub-surface meridional circulation to a depth 
of about 15% of the solar radius (Giles et al. 1997; Braun & Fan 1999). So far there is no 
convincing observational evidence for an equatorward counter-flow deeper down, though 
it must exist to conserve mass. It is generally believed that the turbulent stresses in 
the SCZ drive the meridional circulation, although the details of how this happens are 
not understood (see, for example, Gilman 1986, §3.4.2). So the meridional circulation 
is expected to be confined to the convection zone. However, it is conceivable that an 
equatorward meridional transport of material takes place in the overshoot layer below the 
bottom of SCZ. This view is supported by recent simulations of solar convection (Miesch 
et al. 2000). A dynamo model with a meridional flow through the convectively stable 
overshoot layer seems, at the present time, to be the model with minimal extraneous 
assumptions which gives satisfactory results. 

In the next section, we describe the basic features of our model. Then §3 focuses on the 
parity question and discusses the conditions to be satisfied to ensure the correct parity. In 
§4 we present some more details of what we consider our standard solar dynamo model, 
along with comparisons with observations. Finally our conclusions are summarized in §6. 
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2 Mathematical Formulation 



2.1 The basic equations and boundary conditions 

All our calculations are done with a code for solving the axisymmetric kinematic dynamo 
problem. An axisymmetric magnetic field can be represented in the form 



B = B(r, 9)e (j) + V x [A(r,9)e <t> ], 



(1) 



where B(r, 9) and A(r, 9) respectively correspond to the toroidal and poloidal components. 
The standard equations for the so-called auo dynamo problem are: 



OA 
~dt 



+ -(v.V)(sA) = Vp V 2 - -r )A + aB, 



(2) 



dB 1 

dt r 
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Vt 
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1 dr] t dB 
r dr dr 



(3) 



where s = r sin 9. Here v is the meridional flow, is the internal angular velocity of the 
Sun and a is the coefficient which describes the generation of poloidal field at the solar 
surface from the decay of bipolar sunspots. We allow the turbulent diffusivities r\ v and 
r]t for the poloidal and toroidal components to be different. We describe below how we 
specify v, f2, r] p , rj t , and a. Once these quantities are specified, we can solve (2) and (3) to 
study the behaviour of the dynamo. Apart from the specification of these parameters, we 
also include magnetic buoyancy in a way described below. We carry out our calculations 
in a meridional slab Rb = O.55i? < r < R & , < 9 < it. 

The boundary conditions are as follows. At the poles (9 = 0, n) we have 

A = 0, B = 0. (4) 

For a perfectly conducting solar core, the bottom boundary (r = R b ) condition should be 



A = 0, B = 0. (5) 

However, if the bottom of the integration region is taken well below the depths to which 
the meridional circulation reaches (and hence below the depths to which magnetic fields 
are carried), then the solutions are rather insensitive to the bottom boundary condition. 
We carried out some calculations by changing the bottom boundary condition of the 
toroidal field from B = to 

§~ r (rB) = 0. (6) 
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The solutions remained virtually unchanged. At the top (r = Rq), the toroidal field has 
to be zero (B = 0) and A has to match smoothly to a potential field satisfying the free 
space equation 



Dikpati and Choudhuri (1994) describe how this is done. 

We have used a finer grid resolution as compared to other models existing in the 
literature, with 129 x 129 grid cells in the latitudinal and radial directions. The algorithm 
used by us in developing the numerical code is described in the Appendix of Dikpati 
& Choudhuri (1994) and the Appendix of Choudhuri & Konar (2002). If either the 
a coefficient has a quenching factor or magnetic buoyancy is included to suppress the 
growth of the magnetic field, then any arbitrary initial condition either asymptotically 
goes to zero (sub-critical condition) or relaxes to a steady dynamo solution (super-critical 
condition). We now discuss how the various parameters are specified. 

2.2 Internal rotation Q 

We use the following analytic form to represent the solar internal rotation (Schou et al. 
1998; Charbonneau et al. 1999): 



This analytical expression fits the results of helioseismology fairly closely for the following 
values of parameters: r t = 0.7R & , d t = O.O5i? Qrz/2it = 432.8 nHz, Qscz(@) = ^eq + 
a 2 cos 2 (6) + a 4 cos 4 (6), with n EQ /27r = 460.7 nHz, a 2 /27r = -62.69 nHz and a 4 /2n = 
—67.13 nHz. A contour plot of Q generated by the above expression is shown in Fig. d 
in which the tachocline is shown as a shaded region. Since helioseismology has already 
determined Q, we do not have much freedom to vary its parameters. At the present 
time, however, there exist some uncertainties as to the exact location of the tachocline 
and its thickness. These are indicated by the parameters r t and d t . Especially, it is 
still not completely clear how the tachocline is located with respect to the bottom of the 
convection zone (whether it is completely below the convection zone or partly inside it). 
There are some indications that the tachocline may actually have a prolate shape, such 
that more of the tachocline comes within the convection zone at higher latitudes than 
at the lower latitudes. Given the many other uncertainties in the problem, we have not 
taken this effect into account in our calculations. 
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BCZ (= 0.71 R) 



Figure 1: Contours of constant angular 
velocity as obtained from (8) and used in 
all of our calculations. The tachocline is 
indicated as the shaded region. The bot- 
tom of the convection zone (BCZ) is at 
O.71R.0. The differential rotation being 
symmetric about the equator, only one 
quadrant is shown. 



Figure 2: Streamlines of meridional cir- 
culation obtained by taking penetration 
radius R p = O.61i? . The tachocline is 
shown clS 8b shaded region. Arrows de- 
note the direction of flow. 



We take the bottom of the convection zone at r = O.71-R , which is marked in Fig. 
By bottom of the convection zone, we mean the depth at which the temperature gradient 
changes from being sub-adiabatic below to super-adiabatic above. As we point out later, 
magnetic buoyancy is assumed to be operative only above the bottom of the convection 
zone. Since the strong toroidal field is generated in the tachocline, a tachocline below the 
convection zone would imply a situation where the toroidal field is created at a location 
which is immune to magnetic buoyancy. Hence, the location of the tachocline with respect 
to the base of the convection zone is of considerable importance in our problem. As seen 
in Fig. ^ we take part of the tachocline above the bottom of the convection zone. All our 
calculations are done with this assumed profile of Q. 

In some of our earlier work (Choudhuri et al. 1995; Nandy & Choudhuri 2001; Nandy 
2002), we had taken Q to be independent of 9 and used a radial variation appropriate for 
the equatorial region. If Q is taken to be a function of r alone, then the behaviour of the 
dynamo is much simpler. On the other hand, when Q is a function of both r and 9 as 
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Figure 3: Plot of ^ (in m s x ) as a function of r at the mid- latitude 9 = 45°. 



given by (8), the problem becomes immensely more complicated. 



2.3 Meridional circulation v 

We now describe how the meridional circulation is specified in the northern hemisphere. 
The circulation in the southern hemisphere is simply obtained by a mirror reflection of 
the velocity field across the equator. We get the meridional circulation from the stream 
function ip defined through the equation 



pv = V x [ip(r,9)e^\. 



Assuming a density stratification 



(9) 



(10) 



we take 

ipr sin 9 = ipo( r — Rp) sin 



7T (r - R p ) 
(R Q - Rp) 



{1 _ e _/3irr }{l - e ftr(e - 7r/2) }e~ ((r_r ' o)/r)2 (ir 



with the following values of the parameters: Pi — 1.36 x 10 8 m 1 , p 2 — 1.63 x 10 8 
m" 1 , e = 2.0000001, r = (R Q - # 6 )/4.0, T = 3.47 x 10 8 m, 7 = 0.95, m = 3/2. 
It is the parameter R p which determines the depth to which the meridional circulation 
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penetrates. Fig. El shows the streamlines of the meridional circulation obtained by taking 
R p = 0.61.Rq, which is used in all our calculations and which corresponds to the meridional 
flow going slightly below the tachocline. The amplitude of the meridional circulation 
is fixed by ipo/C. We choose it such that the poleward flow near the surface at mid- 
latitudes peaks typically around Vq = 22 m s _1 . The equatorward counterflow peaks at 
the core-convection zone interface and has a value of 1.8 m s _1 , which is similar to what 
observational analysis of sunspot drift suggests (Hathaway et al. 2003). Since the form of 
the meridional circulation seems crucial for the stability for the dynamo, let us make some 
comments on it. On comparing (11) with equation (9) of Nandy & Choudhuri (2001), it 
will be found that we have now added an extra factor (r — R p ) just after ip . This extra 
factor ensures that vg smoothly falls to zero at R p . If this factor is not included, then vg 
has a finite value at R p , leading to a discontinuity at R p if there is no flow below. Fig. El 
shows the profile of v g at the mid-latitude obtained with the factor (r — Rp) included. It 
should also be noted in Fig. |3] that vg in our model decreases monotonically below the 
bottom of the SCZ and becomes very small in the tachocline. We want to emphasize that 
we need very modest flows below the tachocline (presumably not much beyond the region 
of overshooting) to make our model work. 

We point out that we find well-behaved periodic solutions only for certain forms of the 
meridional circulation. Kiiker et al. (2001) also found oscillatory solutions only within 
limited regions of parameter space. However, the dynamo becomes very robust and stable 
with a sufficiently deeply penetrating meridional flow having a smooth vg profile. 

2.4 Diffusion coefficients rj p and rj t 

We expect the turbulent diffusivity inside the convection zone rjscz to be much larger than 
the diffusivity rjnz in the radiative interior, the overshoot layer being the region within 
which the value of diffusivity makes a transition from r^z to rjscz- Since the poloidal 
component is weak, the turbulent diffusivity acts on it without any difficulty We take 
the diffusivity for the poloidal component to be 



FigE] shows a plot of r/ p with the following values of the parameters: r/scz = 2.4 x 10 12 cm 2 
s _1 j Vrz = 2.2 x 10 8 cm 2 s -1 , tbcz = 0.7R Q , d t = 0.05R Q . The toroidal component, 
however, has a value larger than the equipartition value till it rises to about 30,000 km 
or 40,000 km below the solar surface (Longcope & Choudhuri 2002, §2). Only in the top 
layers of the SCZ, the diffusion coefficient of the toroidal component should be equal to 
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Figure 4: Plots of rj p (r) (solid) and rjt(r) 
(dashed) as given by (fT2|) and (fTT?j) as 
functions of the fractional radial dis- 
tance (r/i? ). y-axis is in units of cm 2 
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Figure 5: Plots of rj p (r) (solid) and rj t (r) 
(dashed) in cm 2 s -1 used in §3.2, as 
function of the fractional radial distance 
(r/R Q ). 



the usual turbulent diffusivity value r\scz- Within the main body of SCZ, the action of 
turbulent diffusivity on the toroidal component must be considerably suppressed. In view 
of this, we take the diffusivity T]t of the toroidal component as shown in Fig. HI which is 
generated from the expression 



r} t {r) 



with rjscz 



Vrz + 



Vsczi 



1 + erf 



BCZ 



+ 



Vscz 



1 + erf 



rTCz\ 



d t ) 



(13) 



4 x 10 10 cm 2 s~\ 



BCZ 



O.72i? and r TC z = 0.95i? Q . As we shall see 



below, rjp and rj t specified in this way gives solutions with dipolar parity. Except §3.2-3, 
everywhere else in our paper we use r\ v and rj t as specified above. 

Simulations of the evolution of the weak, diffuse field on the solar surface (Wang 
et al. 1989a,b), as well as observational estimates from sunspot decay, point out that 
rjscz must be of the order of 10 12 cm 2 s -1 in the upper layers of the convection zone. 
We have taken a value of this order on the higher side and find that it allows sufficient 
diffusion of the poloidal component across the equator to enforce the dipolar mode. A low 
value of diffusivity below the bottom of the convection zone is very important in dynamo 
models with meridional flow penetrating below the tachocline. A low diffusivity in the 
tachocline and the overshoot layer ensures that the toroidal field which is produced in 
the high latitudes within the tachocline does not decay much while being transported to 
low latitudes by the meridional flow (Nandy 2002). The assumed value of tjrz essentially 
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ensures that the magnetic field is frozen for time scales of the order of dynamo period. 
By taking a low r) t at the bottom of SCZ, we also make sure that there is not much 
cross-diffusion of the toroidal component across the equator and it is possible for the 
toroidal component to have opposite values in the two hemispheres. It may be noted that 
there is a term involving drj t /dr in the evolution equation (3) for the toroidal component. 
This term has the form of an advection term, with drj t /dr corresponding to a downward 
velocity. We discovered an error in the original code which produced the results of Nandy 
& Choudhuri (2002). The error made this term involving drjt/dr somewhat smaller than 
what it should have been. However, on incorporating the suppression of rj t in the body of 
the SCZ as we do here, the gradient drjt/dr is moved to the upper layers (which can be 
seen from Fig. QJ and the results essentially remain the same as earlier. 

For the sake of comparison, we present in §3.2-3 some results obtained with a lower 
diffusivity for the poloidal field. The profile of this r\ is shown in Fig. E] This low 
diffusivity does not allow the poloidal components in the two hemispheres to connect 
across the equator and usually the quadrupolar mode is preferred, as we shall see in 
53.2-3. 



2.5 The a coefficient 



The a coefficient is taken in the form 

.1 



a = «o cos9- 



1 + erf 



r — T\ 
d\ 



x 



erf 



r-2 



di 



(14) 



The parameters we use are r\ = O.95_R , r 2 = R Q , di — d% — 0.025i? Q , making sure 
that the a effect is concentrated in the top layer 0.95i? Q < r < R Q . The solid line 
in Fig. El shows the variation of a with r. The amplitude «o is taken such that the 
dynamo is super-critical. We had taken a = 25 m s -1 in most of our calculations, 
where the solutions are found to be super-critical for such a value of «o- We use this 
form of a coefficient in all our calculations except in §3.3, where an additional a effect 
within the SCZ is included. The dashed line in Fig. El shows the radial profile of a used 
in §3.3, the angular variation being taken as cos 9 as used elsewhere in the paper. It 
should be noted that a coefficient in a Babcock-Leighton dynamo is not given by the 
mean helicity of turbulence as in mean field MHD (see, for example, Choudhuri 1998, 
§16.5). In our approach as well in the approach of several other authors (Choudhuri et 
al. 1995; Dikpati & Charbonneau 1999; Nandy & Choudhuri 2001; Kiiker et al. 2001), 
the a coefficient phenomenologically captures the effect of poloidal field generation from 
the decay of tilted active regions near the solar surface. The angular factor cos 9 arises 
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Figure 6: The Babcock-Leighton a used for most of our calculations (solid line) and also 
in our standard model described in §3.1. The dashed line shows the a used in §3.3, as a 
function of fractional radial distance r/R & . Units are in m s _1 . 

from the angular dependence of the Coriolis force which causes the tilts of active regions. 
Several groups suppressed a artifically at high latitudes to reduce the strength of the 
toroidal field at high latitudes. Dikpati & Charbonneau (1999) used the angular factor 
cos 9 sin 9, whereas Kiiker et al. (2001) used cos 9 sin 2 9. Flux tube simulations which 
calculate the tilts of emerging active regions (D'Silva & Choudhuri 1993; Fan et el. 1993) 
give us some clues about the angular dependence of the a coefficient. If the Coriolis 
force could freely produce the tilt without any opposing force operative, then the tilts 
at different latitudes would have been proportional to the Coriolis factor cos 9 and the 
a coefficient would have the same angular dependence. However, the Coriolis force is 
opposed by magnetic tension which tries to prevent the tilt from becoming very large. 
We thus find that the tilt does not increase with latitude as fast as the Coriolis factor 
cos 6*. On these grounds, we expect that a will not increase with latitude as fast as cos 9, 
but taking the angular dependence to be cos 9 sin 9 or cos 9 sin 2 9 may be unrealistic. 

We also point out that we have eliminated the a-quenching factor typically taken to 
be of the form (1 + B 2 /Bq)~ 1 . If magnetic buoyancy is not included, then such a factor 
is the only source of nonlinearity in the model and is essential to allow the simulation 
to relax to steady solutions. We found the a-quenching to be redundant in presence 
of magnetic buoyancy which limits the toroidal field to values < 10 5 G (see §2.5). On 
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theoretical grounds also, the removal of a-quenching is quite logical in the Babcock- 
Leighton dynamo models. In mean field MHD, the a effect comes from helical turbulence 
which is quenched when the magnetic field is super-equipartition. In our model, however, 
the a effect is due to the decay of tilted bipolar regions. Flux tube simulations do show 
that the tilt is less for stronger magnetic fields at the bottom of the SCZ (see D'Silva & 
Choudhuri 1993). So the a effect should depend on the magnetic field at the bottom of 
SCZ, but not on the local value of the magnetic field at the surface where the a effect 
is operative. Since our formulation of magnetic buoyancy (see §2.6) makes the toroidal 
field erupt only when it is of order 10 5 G, we do not expect much variation in a with the 
magnetic field. 

2.6 Magnetic buoyancy 

We prescribe magnetic buoyancy in a way which has been discussed in detail by Nandy 
& Choudhuri (2001) and Nandy (2003). We search for toroidal field B exceeding the 
critical field B c = 10 5 G, above the base of the SCZ taken at r = 0.71 at intervals of 
time t = 8.8 x 10 5 s. Wherever B exceeds B c , a fraction / = 0.5 of it is made to erupt 
to the surface layers, with the toroidal field values adjusted appropriately to ensure flux 
conservation. The parameter / measures the strength of magnetic buoyancy. It was found 
by Nandy & Choudhuri (2001) that, when / was still small compared to 1 (/ has to be less 
than 1), magnetic buoyancy already reached saturation and the character of the dynamo 
did not change any more on increasing /. The value / = 0.5 used throughout our paper 
already puts the dynamo in the buoyancy-saturated regime. 

Although A and B in general evolve according to (2) and (3), we allow abrupt changes 
in B after intervals of r to take account of magnetic buoyancy. While our treatment of 
magnetic buoyancy may not be fully satisfactory, note that uncertainties remain in the 
way buoyancy has been handled earlier by other groups, e.g. by treating buoyancy as 
a simple loss term (Schmitt & Schiissler 1989) or by treating it in a non-local manner 
by making the poloidal source term near the surface proportional to the toroidal field 
strength at the bottom of the SCZ (Dikpati & Charbonneau 1999). 

3 The parity question 

All our calculations are done with internal rotation fl, meridional circulation v and mag- 
netic buoyancy as specified in §2.2, §2.3 and §2.6. We first present results in §3.1 obtained 
with the diffusion coefficients shown in Fig.|3]and an a effect concentrated near the surface 
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as shown in Fig. |H1 by the solid line. The dipolar mode is the clearly preferred mode and 
the results qualitatively match the observational data quite well. Some authors (Dikpati 
& Gilman 2001; Bonanno et al. 2002) obtained anti-solar quadrupolar modes in their cal- 
culations. We believe that this was due to the low diffusivity of the poloidal component 
which did not allow this component to connect across the equator. We present results in 
§3.2 in which the diffusion coefficients are as shown in Fig. i.e. the diffusivity of the 
poloidal component is reduced compared to what we use in §3.1 and that of the toroidal 
field is increased. We end up with the anti-solar quadrupolar parity in this case. Then 
we show in §3.3 that we get back the dipolar parity if we include an a effect within the 
SCZ (as shown by dashed line in Fig. EJ), while keeping all the other things the same as 
in §3.2. This is again in agreement with what Dikpati & Gilman (2001) and Bonanno et 
al. (2002) found. We are thus able to reproduce the results of Dikpati & Gilman (2001) 
and Bonanno et al. (2002). However, we do not agree with their conclusion that a pure 
Babcock-Leighton a effect concentrated near the solar surface cannot give the dipolar 
parity. As we show in §3.1, the dipolar parity is the preferred parity if the diffusivity of 
the poloidal component is sufficiently high, even when the a effect is concentrated near 
the solar surface. 

3.1 Solution with dipolar parity 

We first present a purely Babcock-Leighton dynamo (a effect concentrated near the sur- 
face as shown in Fig. El by solid line), which settles into dipolar parity. As we discussed 
above, the preferred parity depends on the diffusion coefficients. Various experiments 
with the parameter space of the model tell us that there are two important conditions for 
getting the right parity. 

1. The term tjrz representing the molecular diffusivity in the overshoot layer and the 
radiation zone below must be sufficiently small(~ 2.2x 10 8 cm 2 s _1 ) to prevent the 
toroidal field from diffusing across the equator. 

2. The diffusivity of the poloidal field r\ v within the SCZ must be sufficiently large (~ 
2.4 x 10 12 cm 2 s" 1 ) to allow diffusive coupling of the poloidal field between two 
hemispheres. 

Further, we have to avoid a large gradient dr) t /dr at the bottom of the SCZ in order to 
obtain well-behaved solutions. The particular solution we present here is obtained with the 
diffusion coefficients as given in Fig. EJ To ensure that the dipolar parity is the dominant 
mode in the model, we start from a pure quadrupolar initial condition and find that the 
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Figure 7: A time-latitude plot showing the evolution of the toroidal field at the bottom of 
the convection zone (r = 0.71 R®) starting from an initial quadrupolar (anti-solar) mode 
to a pure dipolar (solar-like) mode taking a duration of 3000 years, for the case presented 
in §3.1. The upper panel shows the evolution from 1000 - 2000 years whereas the lower 
panel traces the last thousand years (2000 - 3000 years). Regions shaded in white show 
positive B, whereas black regions denote negative B. 

solution eventually relaxes to a pure dipolar parity. Fig. [7| shows the evolution of the 
initial quadrupolar field into a dipolar field within a duration of 3000 years. Fig. |S] gives a 
snapshot of the relaxed magnetic field configuration, showing that the poloidal field lines 
have connected across the equator, whereas the toroidal field within the tachocline has 
opposite signs on the two sides of the equator. We carried out the simulation for another 
3000 years after the solution relaxed to the dipolar mode, and the solution remained in 
the solar-like parity state during the entire run. We shall provide further details of this 
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Figure 8: A snapshot of (a) the contours of toroidal field B and (b) streamlines of the 
poloidal field, given by contours of constant AsinO, for the case of dipolar parity solution 
presented in §3.1. The solid lines denote positive B or A, whereas dashed lines denote 
negative B or A. 

solution in §4. 

3.2 Solution with quadrupolar parity 

We now present a solution obtained by keeping all the parameters the same as in §3.1, 
except that we change the diffusion coefficients to what is shown in Fig. |3] rather than 
what is shown in Fig. 0] (as used in §3.1). The lower diffusivity of the poloidal field does 
not allow it to connect across the equator and we find that the quadrupolar mode is 
preferred. To make sure that indeed the dominant mode in this case is the quadrupolar 
mode, we started with an initial condition having dipolar parity and found that it relaxed 
to a quadrupolar parity. Fig. |§] presents a theoretical butterfly diagram showing this 
transition process, whereas Fig. El illustrates the magnetic field configuration after the 
dynamo has settled into a quadrupolar mode. We find that the poloidal field has not 
diffused enough to connect across the equator, but has remained separated on the two 
sides of the equator. On the other hand, the toroidal field in the tachocline has the same 
sign on the two sides of the equator and makes up a patch of a common sign across the 
equator. We have made many runs for other values of diffusion coefficients. As we already 
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Figure 9: Time-latitude plot showing the evolution of B at the bottom of the convection 
zone from an initial dipolar state to the preferred quadrupolar state, taking a duration of 
2000 years, for the case presented in §3.2. Only the last 1000 years of model run is shown. 

mentioned, a low diffusivity tjrz below the bottom of the SCZ is an essential requirement 
to obtain the dipolar parity (to ensure that the toroidal field in the tachocline cannot 
diffuse much across the equator). When we take t/rz larger than about 2 x 10 9 cm 2 s" 1 , 
we found that we always got quadrupolar solutions and it was not possible to get dipolar 
solutions even by increasing rj p in the SCZ to facilitate the coupling of the poloidal field 
across the equator. We have shown in Fig. a transition from a dipolar parity to a 
quadrupolar parity in a case in which the quadrupolar parity is preferred. How fast such 
a transition takes place depends on the value of r\Rz- On using rjnz ~ 2 x 10 10 cm 2 s _1 , 
the change-over takes place within just 300 years, whereas decreasing t]rz by two orders of 
magnitude stretches the time scale of transition to the quadrupolar parity to 2000 years. 

3.3 Solution with a effect inside the SCZ 

Finally we consider what influence an additional a effect within the SCZ has on the parity 
of the solution. For this purpose, we change the radial profile of the a effect to what is 
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Figure 10: A snapshot of (a) the contours of the toroidal field and (b) streamlines of the 
poloidal field, for the quadrupolar parity solution presented in §3.2. 

shown in Fig. El by a dashed line rather than the solid line in the same figure, keeping 
the other things the same as in §3.2 (including the diffusion coefficients which are as 
shown in Fig. EJ). We find that the solution relaxes to dipolar parity even if we start from 
quadrupolar parity, as shown in Fig. [TT] A snapshot of the field configuration is presented 
in Fig. This is a case in which the quadrupolar parity would have been preferred if 
the additional a effect inside the SCZ were not present, as we saw in §3.2. However, 
this additional a effect, even though its value inside the SCZ is much smaller than the 
value of the Babcock-Leighton a at the surface, can make the solution dipolar. This is 
in agreement with what has been found by Dikpati & Gilman (2001) and Bonanno et 
al. (2002). One has to keep in mind that the a coefficient has to be multiplied by the 
toroidal field B to provide the source term for the poloidal field. Since B is very large at 
the bottom of the SCZ, even a very small a there can make aB as large as what it is near 
the surface. That is why we find that even a very small a inside the SCZ or at its base 
can affect the nature of the solution so drastically. Presumably, an a effect within the 
SCZ creates some poloidal field there which can diffuse across the equator more efficiently 
than poloidal field created near the surface, thereby enforcing the dipolar parity. Note 
that in the presence of a small a inside the SCZ the magnitude of the Babcock-Leighton 
a required at the surface for steady oscillating solutions decreases by almost a factor of 
10 (dashed line in FigEJ). 
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Figure 11: Time-latitude plot showing the evolution of the toroidal field B at the bottom 
of the convection zone, with an additional a introduced inside the SCZ, as presented in 
§3.3. Last 1000 years of the 2000 year long model run is shown. 

Choudhuri (2003b) has argued that the strong toroidal field at the bottom of the 
SCZ would be highly intermittent and the a effect may be operative in the intervening 
regions of weak field. Several other authors (Ferriz-Mas et al. 1994; Dikpati & Gilman 
2001) argue that the various instabilities associated with the strong field also can produce 
something like an a effect. So an additional a effect within the SCZ or at its bottom 
is certainly a realistic possibility. However, our knowledge about it at the present time 
is very incomplete and uncertain. On the other hand, we see tilted active regions decay 
on the solar surface and we directly know from observations that there is an a effect at 
the solar surface. By the Occam's razor argument, we feel that it is desirable to first 
construct solar dynamo models with this Babcock-Leighton a alone — especially since we 
have shown in §3.1 that it is possible to get solar-like dipolar parity with such dynamo 
models. We present a more detailed discussion of such pure Babcock-Leighton dynamo 
models in §4. 

4 Towards a Standard Model 

We have seen in §3.1 that a CDSD model with an a effect concentrated near the surface 
and with appropriate values of various parameters can give a solution with the dipolar 
parity. We would refer to this solution as our standard model. We have focused primarily 
on the parity of this solution in §3.1. Now we discuss other aspects of this solution and 
show that this solution matches observational data quite well. We have already discussed 
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Figure 12: A snapshots of (a) the toroidal field contours and (b) streamlines of the 
poloidal field, for the solution presented in §3.3. 

in §2 and §3.1 how the various parameters of this particular case are specified. The 
amplitude of the Babcock-Leighton a effect used to generate our standard solution is 25 
m s -1 . For this standard case, we refined our grid to have 257x257 points, and the results 
remained unchanged with this finer resolution. 

Fig. EH shows the time-latitude contour plot of the radial field at the solar surface, 
with the theoretical butterfly diagram superimposed upon it. The butterfly diagram is 
producing by marking the locations of eruption, '+' indicating the positive value of B at 
the bottom of the SCZ which erupts and 'o' indicating the negative value. The sunspot 
eruptions are confined within ± 40° and the butterfly diagrams have shapes similar to 
what is observationally found (see, for example, Fig. 6.2 in Choudhuri 2003a; Hathaway et 
al. 2003). The weak radial field migrates poleward at higher latitudes, in conformity with 
observations. One of the important aspects of observational data is the phase relation 
between the sunspots and the weak diffuse field. See §1 of Choudhuri & Dikpati (1999) 
for a detailed discussion of this (especially note Fig. 1 there). The polar field changes 
from positive to negative at the time of a sunspot maximum corresponding to a negative 
toroidal field B at the base of SCZ. This is clearly seen in the theoretical results shown 
in Fig. ED I n FigE] we show four snapshots of the toroidal field contours and poloidal 
field lines taken successively at an interval of l/8 th the solar cycle period which happens 
to be about 25 years in this case. In all the snapshots we see that the poloidal field 
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lines are symmetric about the equator. The rj p in the convection zone is sufficiently high 
(~ 2.4 x 10 12 cm 2 s _1 ) to ensure that the poloidal fields connect smoothly across the 
equator. On the other hand, the toroidal fields on the two sides of the equator, which 
have opposite signs, cannot diffuse together due to the low r/ t near the base of SCZ. 
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Figure 13: Theoretical butterfly diagram of eruptions for our standard model. The 
background shows contours of diffuse radial field. Eruption latitudes are denoted by 
symbols 'o' and '+', indicating negative and positive toroidal field respectively. The 
dashed contours are for negative B r , whereas the solid contours are for positive B r . Note 
that negative toroidal fields give rise to negative radial field near the poles after decaying 
and vice versa in accordance with Hale's polarity law. 

4.1 The effect of velocity quenching 

One important question is whether the equatorward meridional circulation at the base 
of SCZ should be able to advect the strong toroidal field, working against the magnetic 
tension. It has been argued by Choudhuri (2003b) that this should be possible if the 
strong toroidal field is highly intermittent. However, Choudhuri (2003b) concluded that 
the meridional flow may be barely strong enough to advect the toroidal field. If B becomes 
larger than some critical value, then the meridional flow may not be able to carry it. We 
try to capture this effect by checking at intervals of 10 days if the toroidal field exceeds 
a critical value of 1.5xl0 5 G in a region of thickness O.lli?© below a depth of 0.73i? Q . 
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Whenever B is larger than this critical value, we reduce the velocity at that grid point by 
a factor of eight. A butterfly diagram similar to Fig. EI] is produced in this case and is 
shown in Fig. On comparing Fig. Eland Fig. we find that this velocity quenching 
improves the appearance of the butterfly diagram. We saw in Fig. that eruptions for 
a new half-cycle began at high latitudes before the eruptions at low latitudes stopped for 
the previous half-cycle, leading to a tail-like attachment in the butterfly diagram. We see 
m Fig. CHI that this is gone and a new half-cycle begins at the high latitude at about the 
time when the old half-cycle dies off at the low latitude — in conformity with observational 
data (see Fig. 6.2 of Choudhuri 2003a). 

5 Conclusion 

We show that a CDSD model with an a effect concentrated near the surface and a merid- 
ional circulation penetrating below the tachocline provides a satisfactory explanation for 
various aspects of the solar cycle. The helioseismically determined differential rotation is 
strongest at high latitudes within the tachocline and there is little doubt that the strong 
toroidal field should be produced there, to be advected by the meridional flow to low 
latitudes where it erupts. This view, which was first put forth by Nandy & Choudhuri 
(2002), is a departure from the traditional viewpoint that the toroidal field is produced 
basically at the same latitude where it erupts. Ironically, this new viewpoint makes the 
original motivation of Choudhuri et al. (1995) in introducing the meridional circulation 
somewhat redundant. A standard result of dynamo theory (without any flow in the merid- 
ional plane), which was first derived by Parker (1955), is that the product of a and dQ/dr 
should be negative in the northern hemisphere, to ensure the equatorward propagation 
of the dynamo wave (see, for example, Choudhuri 1998, §16.5). The tilts of bipolar re- 
gions on the solar surface suggest that a should be positive in the northern hemisphere, 
whereas helioseismology found dQ/dr also to be positive in the lower latitudes. If the 
strong toroidal field is produced within the tachocline at low latitudes where it erupts, 
then the simple sign rule would suggest a poleward propagation. Choudhuri et al. (1995) 
introduced the meridional circulation primarily to overcome this tendency of poleward 
propagation, forcing the dynamo wave to propagate equatorward. If the toroidal field 
is produced at high latitudes where dQ/dr is negative, then the dynamo wave should 
propagate poleward even in the absence of meridional circulation. Although the origi- 
nal motivation of Choudhuri et al. (1995) in introducing the meridional circulation may 
no longer be so relevant, it has become increasingly clear in the last few years that the 
meridional circulation plays a crucial role in the solar dynamo. There are indications that 
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the meridional circulation may actually be the time-keeper of the solar cycle (Hathaway 
et al. 2003). We hope that within the next few years helioseismology will discover the 
equatorward return flow of meridional circulation and may even be able to establish if 
this return flow really penetrates below the tachocline, as required by us. 

Dikpati & Gilman (2001) and Bonanno et al. (2002) had earlier argued that a pure 
Babcock-Leighton dynamo with a concentrated near the surface may not give the correct 
dipolar parity. We have clearly demonstrated that this is not the case. If the poloidal 
field has sufficient diffusivity to get coupled across the equator, whereas the toroidal field 
is not able to diffuse across the equator (since turbulent diffusion is suppressed for the 
strong toroidal field), then we find that the dipolar mode is preferred. We saw in §3.3 
that an additional a effect in the interior of SCZ would help in establishing dipolar parity. 
However, in view of the fact that our knowledge about such an a effect is very uncertain, 
we felt that it is first necessary to study pure Babcock-Leighton dynamo models with a 
effect concentrated near the surface alone. Accordingly, we have taken such a dynamo 
model which gives the correct dipolar parity as our standard model. We have shown in 
§4 that this standard model explains many aspects of observational data very well. We 
are right now exploring whether this standard model can explain some other aspects of 
observational data not discussed by us here. For example, active regions in the northern 
hemisphere are known to have a preferred negative helicity. A theoretical explanation 
for this has been provided by Choudhuri (2003b). Our preliminary investigations based 
on the idea of Choudhuri (2003b) show that our standard model presented in this paper 
would give the right type of helicity. We are now carrying out more detailed calculations, 
which will be reported in a forthcoming paper. We may mention that our code can be used 
for other MHD calculations besides the solar dynamo problem. A modified version of the 
code has been used to study the evolution of magnetic fields in neutron stars (Choudhuri 
k Konar 2002; Konar & Choudhuri 2004). 

Acknowledgements. Most of our calculations were carried out on Sankhya, the parallel 
cluster computer in Centre for High Energy Physics, Indian Institute of Science. D. N. 
acknowledges financial support from NASA through SR&T grant NAG5-11873. 
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t = T/8 




Figure 14: Four snapshots of the toroidal field contours (left panel), and the poloidal 
field lines (right panel) separated by l/8 th of the dynamo period T. The case t = is 
shown in Fig. |HJ The line styles are same as in Fig |H1 
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Figure 15: Theoretical butterfly diagram of eruptions for the case presented in §4.1. The 
various conventions used are the same as in FiglT^l 
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